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Abstract 

The supply function equilibrium (SFE) is a model for competition in markets where 
each firm offers a schedule of prices and quantities to face demand uncertainty, and has 
been successfully applied to wholesale electricity markets. However, characterizing the 
SFE is difficult, both analytically and numerically. In this paper, we first present a 
specialized algorithm for capacity constrained asymmetric duopoly markets with affine 
costs. We show that solving the first order conditions (a system of differential equa¬ 
tions) using spline approximations is equivalent to solving a least squares problem, 
which makes the algorithm highly efficient. We also propose using splines as a way to 
improve a recently introduced general algorithm, so that the equilibrium can be found 
more easily and faster with less user intervention. We show asymptotic convergence 
of the approximations to the true equilibria for both algorithms, and illustrate their 
performance with numerical examples. 
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1 Introduction 


With the presence of demand uncertainty, hrms may choose to compete in supply functions 
- a schedule of prices and quantities that correspond to different realizations of the demand. 
Such concept was hrst introduced by Klemperer and Meyer [Tj, and they named the non- 
cooperative Nash Equilibrium of this type of games the Supply Function equilibrium (SEE). 
Soon, people found that the competition in the deregularized wholesale electricity markets 
bear high resemblance to this formulation, and Green and Newbery [2] hrst applied this model 
to the England and Wales market. Since then, modeling behaviors of wholesale electricity 
markets has been an important application of the SEE. 


The SEE model has attracted tremendous attention from both industry and academia. De¬ 
spite its popularity, people found the SEE model difficult due to the following issues: (1) 
The hrst order necessary conditions of the SEE is a system of ordinary diherential equations 
(ODE), shown in Klemperer and Meyer |T|, but when people try to solve this system of 
ODEs, they usually hud that the solutions are not increasing function^j, which is a require¬ 
ment for feasibility; (2) there can be an inhnite number of supply function equilibria, leading 
to an equilibrium selection problem; (3) it is hard to incorporate capacity constraints and 
general cost functions to the framework, and allowing supply functions to have a free form 
makes the problem even more complicated, therefore many studies are limited to symmetric 
hrms and/or restraining the solution space to functions of simple forms, such as lineaiO and 
quadratic functions. 


Despite all the difficulties, researchers have made substantial progress both in theoretical and 
computational analysis of the SEE. Klemperer and Meyer [1] provided foundational analysis 
of the supply function equilibrium, and compares and contrasts the SEE with the equilibria 

^In this paper, the terms “increasing” and “non-decreasing” are used interchangeably. Supply functions 
must be non-decreasing (increasing) to be feasible, but not necessarily strictly increasing. 

^We use “linear” for first-degree polynomials, which are also called affine functions in some literature. 
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of Cournot and Bertrand games. It also showed the existence of the SFE of symmetric 
oligopolies (assuming convex costs, a concave demand, and no capacity constraints), and 
showed that there could usually be inhnite supply function equilibria, unless the support of 
the distribution of the demand shock was unbounded. 

Later, researchers found that capacity constraints could greatly reduce the number of poten¬ 
tial supply function equilibria, and sometimes even make it unique. Holmberg [3] proved that 
if we had symmetric producers, inelastic demand, a price cap and if the capacity constraints 
bound with a positive probability, then we had a unique symmetric equilibrium. With the 
same conditions, except that the producers had identical marginal costs but asymmetric 
capacities, Holmberg jl] showed the equilibrium was unique and piecewise symmetric. 

Perhaps the most important topic in the study of SFEs is how to find them. Finding the SFE 
is difficult, and restrictions are usually needed. Rudkevich et ah [5] provided a closed-form 
formula for cases where demand was inelastic and hrms did not have capacity constraints, and 
it further showed that the equilibrium price had a high mark-up compared to the perfectly 
competitive price. Green [SIT] restricted the supplies to linear functions, and applied the 
model to the England and Wales market. Baldick et ah j8] showed how to hnd linear and 
piecewise linear SFE when the demand and marginal costs were linear. 

A popular approach to Ending the SFE is to work on the first order conditions (a system of 
ODEs). Many of the studies following this approach involved the use of numerical integration, 
but the major difficulty was that the initial conditions were unknown, and without the right 
initial conditions, the integrals so calculated would usually not be increasing functions, i.e. 
they were not feasible supply functions. With the assumption that capacity constraints of 
smaller firms bind earlier, Holmberg |9] provided a procedure for solving the ODE system 
via numerical integration that searched for feasible solutions by tuning the initial conditions, 
i.e., the prices at which the capacity constraints were reached. 
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Baldick and Hogan [TO] proposed an alternative approach that used an iterative scheme 
for hnding the SFE; At each step, each hrm updates its supply function by moving to the 
best response to the other hrms’ previous offers with a discount factor. This procedure 
was repeated and hopefully it would converge. However, as the authors pointed out, the 
computational cost (hnding the best response) of this iterative scheme was huge, even when 
it converged. 

To use the iterative scheme, one must hrst know how to hnd the optimal response to a 
given set of supplies. Anderson and Philpott m provided conditions for the existence of 
the optimal response, and analyzed the bound of difference in proht when approximations 
of the supply functions were needed. Anderson and Philpott [12] and Anderson and Xu 
|13) expressed the expected return of the hrms as line integrals, proved the existence of 
the optimal supply function, and gave necessary and sufficient conditions for optimality. 
Rudkevich da described an algorithm for developing piecewise linear optimal responses by 
cutting the x-y plane into blocks. 

Baldick and Hogan [15] discussed using high degree polynomials in the iterative scheme as a 
parametric form of the supply functions. The authors pointed out that such approximation 
was not stable. By “stable” they meant that given a small perturbation in the equilibrium, 
the supply functions would still converge to the same equilibrium if the hrms followed the 
iterative scheme. 

Anderson and Hu [16] showed conditions for the SFEs’ continuity and diherentiability, which 
served as theoretical guides to algorithm development. In addition, they proposed a numer¬ 
ical method for hnding the SFE that allowed the hrms to have heterogeneous capacities 
and costs. Their method allowed the supply functions to have free form, and approximates 
them with piecewise linear functions. To hnd the equilibrium, the method searched for a 
feasible solution by solving an auxiliary nonlinear program (NLP) that had the necessary 
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conditions as the constraints. This method has been successfully applied in Sioshansi and 
Oren ira. which showed some large generators in the ERGOT electricity market in Texas 
bid approximately in accordance with the SEE model. 

All the above studies of supply function equilibria focused on continuous supply functions, 
while in practice offers in most markets are step-functions. Holmberg et ah |18) showed 
convergence of the discrete SEE to the well studied continuous one as the number of steps 
increases. 

In this paper we beneht from Anderson and Hu [12], and focus on numerical methods for 
hnding the SEE. We allow the supply functions to have free form, and we will exploit the 
capability of splines to approximate the SEE accurately. 

In the hrst part of this paper, we provide a specialized algorithm for markets of asymmetric 
duopolies that have constant marginal costs. In Section |2l we express the hrst order con¬ 
ditions with splines, reduce the problem of solving the system of ODEs to a least squares 
problem, and show that the solution space of the ODE system and that of the least squares 
problem are the same (in terms of approximation). Since we avoided using numerical inte¬ 
gration, we do not have the initial point selection problem. The solution of the least squares 
problem has a very simple form, and we show in Section [3] that we can use a linear search 
to hnd the unique equilibrium of the capacitated problem. In Section 01 we propose an 
improvement to the general method given by Anderson and Hu (TB]. We will see that with 
the use of splines, the number of decision variables and constraints can be greatly reduced, 
thus in principle solving the optimization problem should be easier and faster. Uniform 
convergence will be shown for both methods. Examples that demonstrate the use and the 
properties of these methods are provided in Section [5l 
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2 Solving the First Order Necessary Conditions 


Our model considers a market with m firms. Each firm i has a maximal capacity Capi. Let 
Ci{q) be the cost of firm i for producing an amount of q. Assume Ci{q) is convex, non¬ 


decreasing and differentiable for all i. Each firm knows the exact cost function and capacity 


of its own, as well as those of all the other competitors. 

The market demand is a function of the form D{p, e) = D{p) -t-e. D{p) is strictly decreasing, 
continuously differentiable and concave, and it is known to all firms. The demand shock e: is 
a continuous random variable, and all the firms know that e has positive probability density 
on [SmirD^mai]- We will focus on the type of SEE that each point of the supply function is 
the best response to a realization of £, also termed as “strong SEE” in Anderson and Hu |16) . 
so the knowledge of the exact distribution of £ is not necessary for finding the equilibrium. 

The supply function of firm i is a non-decreasing function Si : [0,pmaa;] —>■ [0,Capj], where 
Vmax = sup{p > 0 I D{p,emax) > 0}. If there is a market specified price cap and if it is less 
than sup{p > 0 | D{p,emax) > 0}, let Pmax equal to the price cap. 

As first pointed out in Klemperer and Meyer [1], in an SEE, the supply functions 
must maximize each firm’s profit 


max p T)(p)+ £ - - Cj{D{p) + e-'^Si{p)), j = 1,... ,m (2.1) 



at all p G [0,Pmaa:]- If fhe supply functions are differentiable at p, and if 0 < Sj{p) < 

Capj for all j, then we have the first order conditions 
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Throughout Section [2] and |3l we assume that the marginal costs are constant. Let the 
marginal cost for firm j be Cj. Then the first order conditions reduce to 



Anderson and Hu [16] proves that in an equilibrium, the supply functions are continuous 
for p ^ Furthermore, it shows that in an equilibrium, each supply function 

Si{p) is continuously differentiable at p ^ {ci,...,Cm} U {pcapi, ■ ■ ■ ,Pcapm}, where pcapi 
is the price where firm i reaches its capacity, i.e., pcapi = inf{p | Sj(p) = Capi}. As¬ 
sume max(ci,..., Cm) < min(pc'api, • • • ,Pcap^)- Since the supply functions are increasing, 
if 0 < Si{p) < Capi is true for all i at a price p, then we must have max(ci,..., c^) < 
p < min(pc-api, • • • ,Pcap,„), which implies that are continuously differentiable at 

p. Therefore, {si\YLi is a solution to the ODE system {\2.2\\ for all the prices p such 
that max(ci,..., Cm) < p < r[im.{pcapn ■ ■ ■ ^Pcap^)- Since we do not know the value of 
min(pcapi, • • ■,Pcap^) yet, we will solve ([22]) numerically for p G {pmimPmax), where pmin = 
max(ci,..., Cm), and we will find min(p( 7 api, ■ ■ ■ ,PCapm) i^i fhe next section. 

Since {sj(p)}™i are continuously differentiable on (max(ci,..., Cm), min(pcapi, • • • ,Pcapm)), 
it is a good idea to approximate them with splines. To achieve continuous differentiability, 
the splines we use should be at least of order 3 (quadratic splines). Order 4 splines (cubic 
splines) are prefered by most people, as they are the lowest-order splines that are smooth to 
human eyes. 

Splines have been very popular for their capability for approximation. And beginning from 
the late 1960’s, splines are being used by mathematicians to develop numerical solutions 
to ordinary and partial differential equations. We will fundamentally do the same in this 
section. To estimate the spline coefficients, one can either use interpolation or use least 
squares estimation. In this paper we use the latter one, and the reason will be justified 
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shortly. 


We start by selecting knots for the spline approximation, and for simplicity, we will let 
the knots for all the supply functions be the same, as this is good enough according to 
our numerical experience. Then, according to the type of splines we use, we will have basis 
functions associated with the knots. Denote the bases with Bt{x), t = 1,..., K, K depending 
on the type and the order of the splines. Let B{x) = {Bi{x ),..., Bk{x)Y', a vector of basis 
functions. Denote the spline approximation of Si{p) with 

Si{p) = '^ktBtip) = B^{p) ■ (3i, 

t 

where bit sue the coefficients to be determined and [3i = (6ji,..., is the coefficient vector 
for Si- 


We replace the supply functions in the first order conditions fl2.2p with their spline approxi¬ 
mations. The equations now become 


'^^huB'tip) 


JZthBtjp) 

p-Cj 


Observe that with p fixed, 02.31) is linear in bit, * 
also be written in matrix form; 


D'ip), j = (2.3) 

1,...,m, t = 1,... ,K. Thus 02.3p can 




( 2 , 4 ) 


where 


= (B'^ip), ..., B'^ip),B'^ip), ..., B'^{p)f, 




p-Cj 


—V— 

m-j 


is a vector of functions, and 












The necessary conditions fl2.4D are linear equations that the splines are expected to satisfy. 
Hence it is natural to use the least squares method to estimate the coefficients, which is 
part of the reason of our choice. At each price p, (12.dh provides m equations. We have 
K ■ m coefficients to estimate, thus one may wish to choose at least K prices from the 

Let the selected prices be 

Pi-)"‘-)PN ^ (.Pmin-) Pmax\ Let 


B = (lB{i}(pi),..., ..., 1 B{i}(p7v), ■ ■ ■, ^{m}iPN)f, 


and 


d = (D'(pi),..., D'(pi),..., D'(p^),..., D'{p^)f 


We expect the spline approximations to satisfy the linear system B/? = d, where a typical 
line of the system, say, ^J^jy{pk)/3 = D'{pk), is a characterization of the relationship between 
Sj and the derivatives of all the other supply functions at price pk- To estimate {bu}, we 
solve the optimization problem 

min ||B/3 — (i||^ . (2.5) 

{bit} 


Before we solve this minimization problem, we would like to have a look at the solution to 
the original ODE system analytically. 


Consider a market with two firms 1 and 2. (12.2p is now a set of two equations: 


■s'lip) 

4(p) 


S2{p) 

P-C2 

Slip) 

p- Cl 


D'ip), 

D'ip), 


^In fact as we will see very soon, it is not enough to determine the coefficients. But to reduce confusion, 
let us just proceed at this point. 
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whose homogeneous problem 


s'i(p) 

4(p) 


S2{P) 

P-C2 

Slip) 

p- Cl 


0 , 

0 , 


has solution 


S2ip) = hip - C 2 ) + 


(c.-c. + (p-c.)log(^)) 

{c2 -Ci + ip- C2) log ^ 


? 


where ti,t 2 ^ that is, a homogeneous solution is a linear combination of two fundamental 
solutions. However, if ^2 7^ 0, it is easy to verify that as p —?■ max(ci, C 2 ) from above, either 
I Si I ^ 00 or |s 2 | —CX3. Therefore for the practical background of our problem, ^2 must be 
0, and consequently the homogeneous solution is just 


si(p) = t(p- Cl), 


S2ip) = t(p- C 2 ). 


( 2 . 6 ) 


Thus if {si(p), S 2 (p)} and {sj(p), S 2 (p)} are two equilibria, we must have sj(p) = s'j’(p) + 
tip —Cl), slip) = S 2 (p) +t(p —C 2 ) for some t. On the other hand, if {si(p), S 2 (p)} is a solution 
to the ODE system (12.2j) . then {si(p) + t(p — ci), S 2 (p) +t(p —C 2 )} is a solution, too, for any 
t G M. Thus fl2.2p has infinite solutions, and we next show that our spline approximation 
can indeed represent all these solutions in duopoly markets. This result holds unless B has 


columns of zeros 


or has fewer rows than columns. 


Theorem 1. For duopoly markets, B does not have full column rank. Furthermore, the 
rank o/B is 2K — 1, unless it contains columns of zeros or has fewer rows than columns. 

^This happens when an knot interval contains no pk if we use B-splines. If we use natural cubic splines, 
it happens if neither of the last two knot intervals contains any pj,. 
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These results do not dependent on the type of splines and the selection of knots. 


Proof. Write the first order conditions in matrix form 




B^{p)/32 

P-C2 


B^{p)/3i 

p- Cl 


+ B^{p)P, 


D'{p), 

D'{p). 


To prove the first part of Theorem [U it is sufficient to show that the matrix of functions 


B^jp) 

\ P-ci 


B^{p) \ 

P-C2 


B'^{p) 




has linearly dependent columns. And since elementary row operations preserve rank, it is 
equivalent to show that the columns of 


/ 

V 


{c 2 -p)B''^{p) B'^{p) 

B'^ip) {ci-p)B'^{p) 


\ 




are linearly dependent. We prove this by using the fact that the elements of B{p) form a 
basis of the space §, which is composed of all the splines on {pmin,Pmax) with the prescribed 
order and knots. 


Assume that we have a nonzero vector = (vf, nj) such that 


/ 


{c 2 -p)B'^{p) B^{p) 


\ 


y B^ip) {ci-p)B^{p) j 


( \ 


Vl 


\V2 ) 


= 0 , 
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or equivalently 


(C 2 - p)B''^{p)vi + B'^{p)v 2 = 0, 

B'^{p)vi + (ci - p)B''^{p)v2 = 0 . 

Let fi{p) = B^{p)vi and f 2 {p) = B'^{p)v 2 - Thus the above can be rewritten as 

{c2-p)f'i{p) + f2{p) = 0 , 

flip) + (ci - p)f 2 ip) = 0. (2.7) 

Recall that flip) and f 2 ip) are splines. Thus (12.7j) implies that flip) and / 2 (p) must be linear 
functions (they are single-piece linear functions because of their smoothness). Therefore, we 
must have flip) = tip — ci) and / 2 (p) = tip — C 2 ), where t is an arbitrary scalar. 

Since flip) = tip — Ci) G § and / 2 (p) = tip — C 2 ) G S, Vi and V 2 must exist and are unique. 
If t 7 ^ 0, then we have Ui 7 ^ 0 and V 2 7 ^ 0, thus we proved that B does not have full rank. 

Furthermore, since t{p — C 2 ) and t{p — Ci) are the only forms that flip) and / 2 (p) can have, 
it implies that the null space of B has only one dimension, i.e., the rank of B is 2K — 1. □ 

Theorem [ 1 ] shows that when we have only two firms, the general solution to the optimization 
problem “minimize ||B/5 — has the form [d = /d*’ -|- f • n, where f G M and where v is 
an eigenvector of B^B whose corresponding eigenvalue is 0. In terms of individual supply 
functions, the general solutions are si = i?^(p)(/3° -f tvi) = B'^{p)(3i + tip — ci) and §2 = 
B'^ip)i/32+tv2) = B'^(p)/32 + t{p — C 2 ), which have the same form as the analytical solutions, 
showing that the splines are able to approximate all the solutions. This is the most important 
reason why we use least squares for the estimation of the coefficients. When the market has 
three or more hrms. Theorem [T] no longer holds — B will generally have full rank, and 
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consequently the optimal solution will be unique and does not represent the solution space 
of the original ODE system. 


We would also like to show the asymptotic property of this spline approximation. We will 
take cubic splines as an illustration. Proofs for other types of splines are essentially the 
same. 


Theorem 2. If Si and $2 are continuously differentiable on \pa,Pb\, where pa > max(ci,C 2 ), 
and if Si and §2 are piecewise cubic splines and are solution to the least squares problem 
h2.5t) . where we place N price levels {pk}k=i uniformly^ among [pa^Vt^ choose 
as the knots, then as the length of the largest knot interval |r| —)• 0 and N —)■ 00 , the ODE 
system will be satisfied by si and § 2 , such that the error functions, s'_fip) — ^^ — D'{p), 
i = 1,2, uniformly converge to 0 on 


Proof. Since Si, §2 and D' are continuous on [pa^Pb], and since Pa > max(ci,C 2 ), the sum of 
squared errors 



Sijp) 

P-Ci 



2 


is Riemann integrable on \pa,Pb\- We will show that the integral 


rPb 

' pa 


E (s-<(p) - - o'iP) 


i=l,2 


P-Ci 


dp —)■ 0, 


as |r| —)■ 0 and W —)■ 00 . 


Since (si,S 2 ) is a solution to the least squares problem, si and S 2 minimize 


N 

E 

k=l 


i=l,2 


(s'-i{Pk) 


HPk) 
Pk Ci 



2 


^This is in fact unnecessary. We place the price levels this way solely for making the Riemann integral 
easier to write. 

^Si is determined by N and the knots {rt}, so it is more rigorous to write Si{p \ N,{Tt}). However, we 
will proceed with Si for conciseness. 
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among all the functions in §, the space that contains all the piecewise cubic splines with the 
selected knots {n}. Let / 4 S 1 and J 4 S 2 denote the complete cubic interpolation of Si and S 2 
with knots {Tt}|3 Since 1451,1452 G §, we must have 


N 

E 

fe=l 

N 

sE 

k=l 


.*= 1,2 


(s'-iiPk) 


.*= 1,2 


SjjPk) 

Pk Q 

hsiiPk 


D'iPk 


Y ( - D\P^ 


Pk Q 


( 2 , 8 ) 


There are upper bounds for the error of the complete cubic interpolations. For example, 
from de Boor [19] one knows that for any g E C^[a, 6], we have \\I 4 g — ^ |r| u!{g\ ■^) 

and \\hg'- g'W^ < ^), where ||/||^ = sup{|/(a;)| | x G [a,b]} and u{f,h) = 

sup{|/(x) — f{y)\ \ \x — y\ < h}. So there exists a constant (F, such that for any g G C^[a, b], 
\\hg - g\\^ < C |r| u{g', and \\hg' - g’\\^ < Cu{g’, ^). 


Therefore as |r| —)• 0, for any p G [Pa,Pb\, 


(hsUp)-b^-D'(p)) 

V p-Ci ) 




5i(p) 


= (/45T(p) - 5T(p)) - 


p-Ci p-Ci 

hsiip) Si(p) 


+ s'_i(p) 


Si(p) 

P-Ci 


D\p) 




p-Ci p-Ci 
2 


—)■ 0 , 


(2.9) 


where the second equality is because s'_^{p) — = D'(p), i = 1,2, a necessary condition 

for an SFE, and the inequality is by applying the error bounds and the triangle inequality. 
Convergence is due to uniform continuity of 5 ', i = 1,2, on [Pa,Pb\- 
^Similar to Si, it is actually l 4 Si{p \ {ti}), but for conciseness we will use /is^. 
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So by using the definition of Riemann integral, we have 


r-Pb 


0 < lim 


U=l,2 


Y.(iUp)-^-D'(p) 


p-Ci 


dp 


= lim lim 

|r|—>0 -/V—>-oo 


< lim lim 

|r|—S'O W —>-00 


Pb-Pa 

N 

Pb-Pa 

N 


N 

E 

k=l 

N 

E 


4 = 1,2 




k{pk 


Pk Q 


D'iPk) 


k=l U=l,2 


( hs'_i{Pk) - 


hsiipk) 


Pk Q 


- D'{p,) 


< lim lim -——N 

|r|—>-0-/V—>-oo N 


4 = 1,2 


—) + 


^ 1 ^ I 
2 


{Pa - Ci) 


= lim {Pb-Pa) 

|r|^0 


4 = 1,2 


" ' ^(^-4, y) + 


7-U , 


I hi' 


{Pa - Ci) 


= 0 , 


where the second inequality is by fl 2 . 8 l) . and the third inequality and the convergence are by 

(ESI). 

The uniform convergence follows naturally as [pa^Pfe] is closed and the error functions s'_^{p) — 
— D'{p), i = 1,2, are continuous on \pa,Pb\- 

When Si and S 2 satisfy stronger conditions, we can use tighter bounds. For instance. Hall |2n) 
and Hall and Meyer [21] show that if S' G C"^[a, b], then the tightest bounds are || J 4 S' — slloo — 
and \\hg' - sIL < ^ □ 


Now we have a simple form of the spline approximations, which we know will converge to 
the true solutions as the mesh of the splines becomes finer. In the next section we will take 
this advantage and find the SFE with capacity constraints. 
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3 SFE of Duopolies with Capacity Constraints 


In Section [2] we solved the necessary conditions for dnopoly markets. In this section, we still 
focus on duopoly markets, and we use the solutions from Section [2] to find the SFE with 
capacity constraints. SFE of more firms will be discussed in the next section. 

In the following, we denote by f{x~) and f{x~^) the left and right limits of /, respectively. 
Similarly, we denote by f'{x~) and f'{x~^) the left and right derivatives of /, respectively. 
Same as in Section |2l we use pcapi for the price where Si reaches the capacity, i.e., pcapi = 
inf{p I Sj(p) = Capi}. 

Proposition 1 . In a 2-firm-SFE, assume firm 1 reaches the capacity earlier than firm 2 
does, i.e., pcapi < Pcap 2 - assume that S 2 {pcapi) > 0. Then si is differentiable atpcapi, 
and the derivative is 0. In other words, the supply function that reaches the capacity first 
must reach it smoothly. 

Proof. Since and S 2 {pcapi) > 0, we have max(ci,C 2 ) < Pcapi < Pcap 2 (see Anderson and 
Hu |I6]). So there exists 5 > 0, such that sfip) is differentiable for p G {pcapi — h,pcapi + 
5)\{pcapi}- And when sfip) is differentiable, the first order condition (12.21) for S 2 {p) can be 
written as 

S2{p) = {p- C2){s[{p) - D'{p)). (3.1) 

Once Si reaches Capi, it cannot decrease, as we require supply functions to be non-decreasing. 
Thus Si(p) = Capi for p > pcapi, and s'fipf^pj = 0. If Si were not smooth at Pcapi, he. 
^liPcapi) > 0, then from (13.1|) . we must have 

S 2 {Pcap^) = (PCapi “ C 2 ) (s'l (p^^^ J - Dfpcapfi) 

> (PCapi - C 2 ){s[{p^^pJ - Dfpcapfi) = S2{Pcap,)- 
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Thus S2 would be decreasing at Pcapi, and it would be disqualified as a supply function. 
Therefore, we must have = s'i(p^^p^) =0. □ 

Proposition 2. Under the same assumptions of Proposition{^ if D[p) is twice differentiable, 
then the derivative of S2 has a jump atpcapi, he., 

Proof. Differentiate both sides of (13.ip . we have 

4 (P) = ( 5 'i(p) - D'ip)) + (p - C2 )(-D"(p) + sfip)). 

Since s'Kpf^pJ < 0, = 0, and Proposition [1] shows that s[{pcapi) = 0, the left and 

right limits of S2 must have the relationship 

4(PCapi) = -Dfpcap,) + (PCap, “ C 2 ) (s'/fepi) “ ^^"(PCapi)) 

< -D'(pcapi) + (PCapi - C2)(s'/(PCapi) “ ^^"(PCapi)) = 4(PCapi)- 


□ 

Propositions [U and [2] help us understand the nature of the SFE with capacity constraints. 
Proposition [1] also provides a hint on how to find the equilibrium. 

Recall from Section |2] that a general solution can be written as ft = fi^ + t-v, where f G M and 
where v is the eigenvector of that corresponds to the eigenvalue 0. Note that /3 is just a 
solution to the ODE system, and B"’"{p) • /5j, z = 1, 2 may well be decreasing or even negative 
at some part of {pmim Pmax) ■ Our aim is to find the by adjusting t, such that B'^{p) -/Si is a 
nondecreasing curve, and has maximum Capi at a price p > Pmin, which we define as pcapi, 
and that B^{p) ■ (32 is nondecreasing from Pmin to pcapi with B'^{pcapi) ■ 1^2 < Cap 2 . (Swap 
1 and 2 if necessary.) If pcapi > Pmax, then the capacities are not binding. If pcapi < Pmax, 
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then Capi is binding, and the estimated equilibrium will be 




max(-L)'(p)(p-ci),0), p < Pmin] 


h{p) = B^{p)-^i, 


Pmin P ^ PCapi i 


Capi, 


PCapi P Pmaxi 


and 




max(-D'(p)(j9 - C2),0), P<Pmin] 


h{p) = B'^{p) • /32, 


Pmin ^ P — PCapi ; 


mm{-D\p){p-C2),Cap2), Pcapi <P< P] 


max • 


\ 


Monotonicity is not an issue at the lower end where p < Pmin = niax(ci,C 2 ). In this price 
range, the hrni with the higher marginal cost does not produce, and the one with the lower 
marginal cost outputs at the monopolistic level. When p reaches max(ci,C 2 ), the high cost 
hrm begins to produce, and fl3.ip tells us that the low cost hrm can only have a sudden 
increase in supply at that price. Hence there is no issue with monotonicity. 

In practice, hnding the appropriate (3 = (3^ + t ■ v is easy. By plotting the splines, we can 
easily spot the trend of how the supply functions change when we adjust t, and we can 
also see intuitively that in general there can be at most one SFE with capacities constraints 
(sometimes an SFE just does not exist). If we are convinced that an SFE exists, we just 
need to do a linear search (thanks to the 1-dimensional solution space) to hnd the t that 
makes one of the supply curves reaches its capacity smoothly, according to Proposition [1] 

Theorem [2] guarantees that in the limit situation, si and S 2 are solution to the ODE system 
fl2.3p for p G {Pmin,Pcapi), and by Proposition 3 in Holmber et ah jlB], the Si and S 2 so 
constructed are indeed a supply function equilibrium. 

Proposition 3. If D{pmin) + ^min < 0, then there can be at most one (strong) SFE with 
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capacity constraints. 


Proof. The condition D{pmin) + ^min < 0 means that it is possible that the demand is 
sometimes really low, and the market clearing price must be lower than Pmin = niax(ci, C 2 ), 
thus by its definition, an SFE has to include prices below Pmin, which further means that when 
p G {pmin, min(p( 7 ap^,P( 7 ap 2 )), the difference between two equilibria has to be t{p—ci) for Si{p) 
and f(p —C 2 ) for S 2 {p), for some f G M, according to fl2.6p . Without loss of generality, assume 
si{p) reaches its capacity first, at pcapi- Proposition [U shows that s[{pcapi) = 0. If si{p) is a 
supply function of firm 1 in any equilibrium, we must have Si{p) = Si{p) + t{p — Ci) for some t. 
If t was positive, then si{p) would reach its capacity at a price p < pcapi- Since s'^(p) > 0, we 
have s[{p) = s[{p) +1 > 0, which contradicts with Proposition [1] If t was negative, then we 
have si(p) = Si{p)+t{p-ci) < Capi for p < pcapi, and s[{pcapi) = s[{pcapi) + t = 0 + t <0, 
which disqualifies si{p) as a supply function. Therefore t has to be 0, which means si = si, 
and hence we cannot have two distinct SFE. □ 

All we discuss in this paper are strong SFEs, which are not guaranteed to exist. However, 
Anderson m shows that at least for duopoly markets, weak SFEs always exist, which is 
beyond the discussion of this paper. 


4 A General Method for Finding SFE 

When the market has more than two firms, the solution to the least squares problem will be 
unique. Thus the method used in Section [3] for finding the SFE with capacity constraints 
will not work, and hence we need a new method. Also, we would like a method that handles 
general cost functions, instead of just linear ones. But first of all, we would like to show 
some conditions that an SFE must satisfy when we have more than two firms. 
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4.1 Properties at the Nonsmooth Points in Multiplayer SFE 


In Section [3] we saw that when we have two firms, Proposition [T] shows that the supply 
function that reaches its capacity first must reach it smoothly. When we have m firms, 
m > 2 , for the same reason, the m-lth supply function to reach its capacity should still 
reach it smoothly, but the first m — 2 supply functions do not have to. 

In Anderson and Hu jlG] the authors show that in an SFE, a supply curve can be discontinu¬ 
ous only at a price where another firm begins to produce, while all the other producing firms 
are at their capacities. This leads to the following consequences, which should be observed 
in a good SFE approximation. 

Suppose Pcapi is the price where Si reaches its capacity Capi, and suppose that firms 1... n, 
n > 2, are producing at Pc^pi ^^t bound by their capacities. Then due to continuity, 

the following must hold; 


1- Siipcap^) = SiiPcapJ = (PCapi “ ^’j^PCapJ “ ^'(PCapJ), * = 1, • • • , U; 

2. Siipcap^) = Si(Pcapi) = (PCapi “ Ci) “ ^'(PCapi)), i + 1- 


And together they imply 


AiVcap^ = (SifepJ - s'(Pcapi)) = ■ ■ ■ = (s'fepi) “ s'(Pcapi)) , 


i^l,2 

which further implies 

- 4(Pcapi) 


27 ^ 1 ,n 



gl(Pcapi) 
n — 2 


It means that the rest of the curves Sj, i 7 ^ 1 are not differentiable at pcap^ and the right 
limits of their derivatives minus the left limits are all equal. Graphically, in an SFE, we 
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expect to see all these curves have a jump in their derivatives at this price by the same 
amount. 

At the lower price level, where hrms begin to produce, similar things happen, but only that 
it is now a decrease in the derivatives: At ci, hrm 1 begins to produce. And suppose that 
hrms 1... n, n > 2, are producing at cj^and are not bound by their capacities. Then we 
must have: 


1 . Si(ci) = Si(c^) = (ci - i = l,...,n; 

2. Si(ci) = Si(ch) = {ci- - D'{c{)), i 7 ^ 1. 


Together they imply 


1 ^ 1,2 




27^1,n 


which further implies 


4(c+)-4(cr) = --- = <(c+)-<(cr) = 


.fM) 

n — 2 


In a graph of the SFE, the already producing hrms will have a drop in their derivatives by 
the same amount, whenever there is a new hrm begins production. 


4.2 A General Method 

In this subsection we develop a general method that works for markets with arbitrary number 
of players, and the cost functions are no longer assumed to be linear. Of course this general 
method can work with duopolies with linear cost functions, but still the method introduced 
in Sections [2] and [3] are recommended, as least squares problems are extremely easy to solve. 
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In Anderson and Hu [16] the anthers show how to use piecewise linear functions to approx¬ 
imate the supply functions in an equilibrium. They list the necessary conditions that the 
supply functions of an SFE must satisfy, and try to find a set of piecewise linear functions 
that satisfy these conditions at selected prices. To do so, they form an auxiliary optimiza¬ 
tion problem with the necessary conditions as constraints, and solve for a feasible solution. 
However, as they report in the paper, a feasible solution is not easy to find. They need 
to relax the equality and inequality constraints to the error being less than or equal to a 
bound, and let the bound shrink to zero with iteration. In addition, this method requires 
user intervention: sophisticated artificial constraints need to be added to help the solver find 
a feasible solution, and according to the authors, some solvers were sensitive to the objective 
function, i.e., under the same constraints, the solver may deem a problem infeasible with 
one objective function, but could find the optimal solution when given an another objective 
function. So when the problem doesn’t solve, the user doesn’t know whether it is because 
the SFE doesn’t exist or it is because he/she is not using the right objective function. 

Here we base on the same idea and improve by simplifying their method with the use of 
splines. In fact, their piecewise linear functions could be seen as splines with free knots (the 
knots were decision variables in their model), but with formal use of splines we can greatly 
reduce the number of variables and constraints of the problem, which in principle makes it 
easier to find a feasible solution, and faster to find the optimal solution. 

If {si}^i form an SFE, then for any firm i, and for any demand shock e, the corresponding 
market clearing price p must solve the optimization problem 

max Profit{p) = [D{p) + e - Sj{p)]p - Ci{D{p) + e - Sj{p)) 
s.t. 0<D{p)+e- E Sj{p) < Capi- 

If are differentiable at the optimal price p, then they must satisfy the following 
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Karush-Kuhn-Tucker (KKT) conditions: 


Si{p) + {p- C'i{si{p)) -Xi + Pi) ■ {D'{p) - T.j^/j{p)) = 0, 

T.jSj{p) = D{p) + e, 

Si{p) < Capi, 

Si{p) > 0 , 

Xi{Capi - Si(p)) = 0, 

PiSi{p) = 0 , 

Ai,/ii>0, (4.1) 

where A* and pi are Lagrangian multipliers corresponding to the capacity and non-negativity 
constraints, respectively. Replace with their spline approximations in (14.Ij) . 

and assemble iu for all firms and a set of demand realizations 5^, k = 1,..., A^. Assume 
that {si\YLi are differentiable at the optimal prices pk = p{£k), then the KKT conditions 
become: 


HPk) + (pk - C-{si{pk)) - Xik + Pik){D'{pk) - Y.j^iS'j{Pk)) = 0, for all i, k, 
Y.jSi{Pk) = D{pk) + Sk, for all k, 

SiiPmax) < Capi, for all i, 

SiiPmin) >0, for all i, 

XikiCapi - hipk)) = 0, for all i, fc, 

Pikhipk) = 0, for all k, 

Xik, Pik > 0, for all i, k, (4.2) 

where Si{pk) is expressed as Si{pk) = B'^{pk)/3i in computation. 

Under what conditions do /3i , Xik and pik exist that satisfy (14.21) ? The answer depends on 


23 





what splines we use and how we select knots. For example, if we use splines of order 4, and 
if we set one knot at each price level, i.e., = p^, then we are sure there exist /dj , 

and pik that satisfy (14.2p . given the SFE itself exists. In fact, conditions (14.111 are all about 
and their first order derivatives. If we have B^{pk)(3i = Si{pk) and B'^{pk)(5i = s'(pfc) 
for all i and k, then the f3i and the original Xik and pik will automatically satisfy fl4.2p . 
This is not difficult. Since = Pk, for any i, from pk to Pk+i, B'^{p)Pi is a single piece 
cubic polynomial. And B'^{pk)(3i = Si{pk), S'^(pfc+i)A = B''^{Pk)/3i = s[{pk) and 

B''^(pk+i)/3i = s[{pk+i) place 4 constraints that will determine the polynomial. A piecewise 
cubic Hermite interpolation is a spline that satisfies these constraints. (See de Boor |19).I 

More generally, if we are using B-splines, for example, then a iF-knot B-spline of order M 
has K + M coefficients. If we consider a price range where 0 < s* < Capi for all i, thus 
all supply functions are smooth and all Xik and pik are 0, then there is only one equality 
constraint per firm per price level, i.e., the first constraint in fl4.2p . Therefore, if we place 
only one price level between every two adjacent knots, then there exist that satisfy 

fl4.2p . However, if we have more than one price level between some adjacent knots, then a 
solution that satisfies (14. 2 p may not exist. 

As in Anderson and Hu [1^, we use an auxiliary optimization problem to find a feasible 
solution. The problem here is that B^{p) has no simple analytical expression, thus it can 
hardly be evaluated by a solver. Fortunately, since {sj} are optimal for all the values of s, 
{sk} do not have to be chosen to reflect the distribution of e. Hence, instead of selecting 
{£fc} and optimizing {pfc}, we can fix {pk}, and let {sk} be the decision variables. Further, 
examining fl4.2D closely, one would find that do not have to appear as decision variables 
at all; they are simply determined by Sk = B'^{pk) /^j ~ D{pk). If there is an Sk larger 

(smaller) than the upper (lower) bound of the support of e, it means that we have chosen a 
Pk too large (too small) that is not needed for the supply functions. 
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In addition, feasible supply functions must be non-decreasing. Thus we place the monotonic¬ 
ity constraint Si{pk) < Sj(pfc_|_i), for all i and k. However, with this new constraint, we are no 
longer guaranteed to hnd /?* , \ik and piu that satisfy (I4.2jl . i.e., there may not be a feasible 
solution in the spline space, which means that we need to do relaxations. We replace the 
“= 0” constraints in fl4.2p with their absolute values less than or equal to p, where p > 0. 

The objective is simply “minimize p”, thus there is no need to use iteration as in |16] . 

The complete formulation is now as follows: 


min p 

P,/3,X,ii 


S.t. 


Si{Pk) + (Pfc - C'{si{pk)) - Xik + Pik){D'{pk) - JZj^i^'jiPk)) 
Si(pfc) < Si{Pk+i), for all i, k, 

SiiPmax) < Capi, for all i, 

SiiPmin) >0, for all i, 

XikiCapi - Si{pk)) < p, for all i, k, 

PikSi{pk) < P, for all i, k, 

Xik,Pik > 0, for all z, k, 


< p, for all z, k, 


(4.3) 


where s,(pr) ^ 

Ideally, we would hope that the optimal value of p to be 0. But in reality, it is often the 
case that there is not a solution in the spline space that exactly hts all the conditions in 
fl4.2p . thus the optimal p will be positive. However, due to the flexibility of splines, there 
are functions in the spline space that “almost” £t fl4.2p . i.e., the optimal p will be small (See 
examples in Section [S]). 

We now show an asymptotic property of the approximation. Let the knots be Pmin = < 
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Ti < • • • < Tat = Pmax, and for simplicity, let the controlled prices be pk = \ijk-i + T'fc), 
A; = 1,..., iV. We will only show and prove the property for quadratic splines, but it can be 
proved very similarly for splines of higher orders. 


Theorem 3. Assume that {si}^i form an equilibrium and are differentiable at {rfc}^ 




and at {pk}k=i that are defined as aboveu Let {si}^^ be the optimal solution to 
among quadratic splines, where the knots are {r^j^Q and prices are {pk}k=i> then every 


price p G [PminjPmax] wUl eventually satisfy the KKT conditions as |r| —)■ 0. 


Proof. We first show that the KKT conditions will eventually be satisfied at all the controlled 
points, i.e., the optimal value p —)■ 0, as |r| —)■ 0. Then we show that the KKT conditions 
will be satisfied at any point p G [pmin^Pmax]- 

For every i, let Sj be a smooth approximation of Sj, such that if [Tk-i,Tk] does not contain 
a nonsmooth point of s,, then Si{p) = Si{p), for all p G [r^-i, r^]; otherwise, we only require 
non-decreasingness and Si{p) = Si{p) at pk- Note that as |r| —)■ 0, we will not have two 
adjacent intervals that both contain nonsmooth points, and each interval will contain at 
most one nonsmooth point. 

Let I^Si be the quadratic interpolation of Sj, such that I^Si{pk) = Si{pk) = Si{pk) for all k = 
1,... ,N. Thus the monotonicity constraint, the capacity constraint and the nonnegativity 
constraint are automatically satisfied. 

Since s* is smooth, we can use the property proved in Marsden [23] that for all k, 

Ihs'iiPk) -s'i{pk)\ < 3sup||s'(p) -s'(p')|: such that \p-p'\ < yj = 0(|r|). 


Let pik and Xik be the same as they are in fl4.ip . so that the complementarity constraints 

®If not differentiable, choose slightly different and {pk}k=i- Remember that every Si{p) is con¬ 

tinuously differentiable only exeept at p G {(71(0),..., (7{,j(0)} U {pcapi, ■ ■ ■ ,PCapm}- 
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are satisfied, and only the first constraint in fl4.3|) will affect the optimal value p. 


Consider [rk-i, r^] that does not contain a nonsmooth point of s*. We have Si{p) = Si{p) for 
all p 6 [Tk-i,Tk]. The error of the first order condition given by I^Si will be: 

hSiiPk) + {Pk - C[{hSi{pk)) - Xik + Pik)iD'{pk) - T.j^ihSj{Pk)) 

= Si{pk) + {pk - C[{si{pk))) - Xik + Pik){D'{pk) - Y.j^ihs'j{pk)) 

= Si{pk) + {pk - Cl{si{pk))) - Xik + Pik){D’{pk) - T^j^iSjiPk) + T.j^i{s'j{pk) - hs'jiPk))) 

= {Pk - Cl{si{pk))) - Xik + Pik){T.j^iis'jiPk) - hs'jiPk))) 

<{Pk - C'{Si{pk))) - Xik + Pik){T.j^i \s'j{Pk) - l3s'j{Pk)\) 

={Pk - Cl{si{pk))) - Xik + Pik){T.j^i\s'j{pk) - hs'j{pk)\) 

<iPk - C'iisiipk))) - Xik + Pik)iJ2j^iOi\r\)) 

=0{\t\), 

where in the first equality we replaced hsiipk) with Si{pk) as they are equal, in the second 
equality we added and subtracted J2jj!=i^'ji.Pk), in the third equality we removed Si{p) + {p — 
C'i{si{p)) — Xi + Pi) ■ {D'{p) — Yhj^i 'Sj(p)) ns it equals 0, and in the fourth equality we replaced 
s'jipk) with s'j{pk), because by construction Sj = Sj on[rfc_i,rfc]. 

Now consider [Tk-i,Tk] that does contain a nonsmooth point p* of Sj. There are 2 cases: 
p* < Pk and p* > Pk (We assumed pk is a differentiable point). 

For the case p* < pk, in the next interval [rk, Tfc+i], for alH = 1,. .., m, (a) s, is smooth, (b) 
Si = Si, and (c) I^Si is a quadratic interpolation of Sj. So we have (d) s){pk) = s[{Tk) + 0(|r|) 
(by a), (e) hs'iijk) = s'(rfc) + 0(|r|) = s'(rfc) + 0(|r|) (by b and c), and (f) hs'iipk) = 
hs'i{Tk) + 0(|r|) (by c). Thus by d, e and f, \s){pk) - hs'i{pk) \ = 0 {t), for all i. 

Similarly, in the case p* > pk, we look at the previous interval. So for all i, s[{pk) = s'(rfc_i) + 
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0(|t|), J 3 s'(rfc_i) = s'(rfc_i) + 0(|r|) = s'(rfc_i) + 0(|r|), and hs'iiPk) = hs’iiTk-i) + 0(|r|), 
thus \s'i{pk) - hSiiPk)\ = 0 {t). 

Therefore, the error of the first order condition given by I^Si will be: 

hSiiPk) + {Pk - C-{si{pk)) - Xik + Pik)iD’{pk) - Y.j^ihs'j{,Pk)) 

<{Pk - C-{si{pk))) - Xik + Pik){T.j^i\s'j{pk) - hs'j{pk)\) 

<iPk - Ci{si{pk))) - Xik + Pik)iJ2jjLiOi\r\)) 

=0{\r\). 


Thus we can conclude that the best p with {/aSj}™! is 0(|r|). And since {/sSi}™! is just 
one of the feasible approximations, the optimal value p for fl4.3p must be at least as good as 
0{\t\). Therefore p —)■ 0, as |r| 0. 

The uniform convergence follows naturally due to uniform continuity: Let Xi{p) and Pi{p), 
i = 1,... ,m, he linear interpolations of {Ajfc} and {pik} (so they are non-negative), i.e.. 




A,; 


iki 


p = Pk] 


Pk-p 

Pk—Pk-i 


Ai,fc-i + 


p-pk-i 

Pk-Pk-i 


X 


iki 


Pk-i <p<Pk, 


and 


Pi{p) = 


f^ik'! P Pki 

Pik-i+ Pik, Pk-i < p < Pk- 

Pk—pk-l^' Pk—Pk-l'^ ' 1 t'K 


So the approximated supply functions {si(p)}, and the Lagrangians {Ai(p)} and {pi{p)} are 
uniformly continuous on \pmin-iPmax\- Also, the error functions of the first order conditions 


FOCi{p) = Si(p) + (p- C'(si(p)) - Xi(p) + pi(p))(B'(p) -Y^Sj(p)), i = 

j¥=i 
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and the error functions of the complementarity conditions 


LAMi{p) = Xi{p){Capi - hip)), i = 1,.. .,m 


and 

MUi{p) = fii{p)k{p), i = 

are uniform continuous. 

For any p G [Pmin,Pmax], let Pk be the nearest point to p among {pk\- Hence, for any e > 0, as 
proved above, there exists > 0, such that when |r| < (5i, for all i, we have \FOCi{pk)\ < f) 
\LAMi{pi)\ < I and \MUi{pi)\ < f- Convergence at all {pk} are dominated by p, so 
is independent of Pk- By construction, 0 < Si{pk) < Capi. Also, there exists 82 > 0, such 
that when \p — p\ < 82 , due to uniform continuity, we have \FOCi{p) — FOCi{p)\ < |, 
\LAMi{p) - LAMi{p)\ < |, \MUi{p) - MUi{p)\ < | and \hip) - Si{p)\ < e. Thus, for 
|r| < min(5i,52), we must have \FOCi{p)\ < e, \LAMi{p)\ < e, \MUi{p)\ < e and —e < 
Si{Pk) ^ Capi + e. Therefore, the KKT conditions will eventually be satished uniformly at 
all p G \pmin,Pmax], aS \t\ 0. □ 

The hrst constraint in (14.31) from the KKT conditions is just what the ODE system f|2.2l) 
says, thus Proposition 3 in |1H] and Theorem [3] together guarantee that in the limit situation 
the solution of fl4.3p is a supply function equilibrium (when it exists). 

In case the solution gets stuck at a local minimum, one can try replacing the constraint 
Siipk) < Si{pk+i) with bi^t < if one uses B-splines. The constraint bi^t < bi,t+i is a 

sufficient condition for non-decreasingness for B-splines, so using it instead of the necessary 
condition Si{pk) < Si{pk+i) will reduce the space of feasible solutions, thus making the 
solution less likely to fall into local minima. Also, one does not want to over reduce the 
space, so quadratic splines are recommended, because bi^t < is a necessary and sufficient 
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condition for quadratic splines to be non-decreasing. In addition, our experience shows that 
although the optimal p given by the pointwise monotonicity constraint Si{pk) < Si{pk+i) is 
usually smaller than that from the full monotonicity constraint bi^t < bi,t+i, the solution from 
the latter formulation is usually more robust than the former (See Example [3]). Thus it is 
always good to consider using the full monotonicity constraint, even when local minimum is 
not present. 

We see that compared to the formulation given in Anderson and Hu (see [1^ (14), (17) and 
(18)), formulation (14. 3 h has significantly less variables and constraints, so in principle it is 
easier to hnd a feasible solution and faster to solve. It is also user friendlier, as it does not 
require the user to adjust the objective function and constraints when solving a problem. 
Since there is no description in [TB] about in which cases it is hard to hnd feasible solutions, 
we are unable to make a comparison, but we have not experienced any difficulty in the many 
problems we tested. 

When solving the problem we selected hnite {pk}- If the hrms were able to choose p contin¬ 
uously, they may be able to improve the proht slightly. The improvement of hrm Ts proht 
at Pk can be approximated by \Profit'{pk) ■ {p — Pk) \ + 0{\Tf). If Xik = 0 and pik = 0, then 
\Profit'{pk) ■ {p — Pk)\ < P I'^l = 0{\Tf ). So the improvement shrinks to 0 as |r| —)• 0. 


5 Numerical Examples 

The following are a few examples demonstrating the use of the numerical methods we intro¬ 
duced above. Without loss of generality, constant terms of all the cost functions are set to 
0, as they do not affect the results. When solving problems with the general method, both 
IPOPT 121] and CONOPT [25] are good choices for the solver. With their default settings, 
CONOPT tends to give a slightly better solution in terms of optimality and feasibility, while 
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IPOPT is much faster. 


Example 1. In this example we use the least squares method to hud the equilibrium in 
a duopoly market. The two hrms have linear cost functions: C'i(g) = lOg, C 2 {q) = 15g. 
Their capacities are Capi = 80 and Cap 2 = 75, respectively. The demand function is 
D{p) = —3p + £. 

We use natural cubic splines in the example, while B-splines work hne, too. The knots are 
from 5 to 77 at step 9. The price levels used for htting the hrst order condition (12.2j) are 
from 16 to 65 at step 0.5. As described by Theorem [H the B matrix has 18 columns but the 
rank is 17, giving us one degree of freedom, which covers all the potential solutions when 
Emin is low enough. Figure 15. lal plots the solutions to fl2.2p with different values of t. 

From Figure IS.lal it is easy to tell that Firm 1 will reach the capacity hrst. A linear search 
gives that the maximum of Si will equal to Capi (at Pcapi = 31.65) when t = 194.06. 
Therefore, the obtained splines with t = 194.06 gives an approximation of the SFE for 
15 < p < Pcapi- When p > Pcapn Si(p) = 80 and S 2 {p) = 3(p — 15) until S 2 reaches Cap 2 
at p = 40. When 10 < p < 15, si(p) = 3(p — 10) and S 2 (p) = 0. Figure [5Tb] shows the 
approximated SFE for 10 < p < 45. 


Example 2. This time we hnd the SFE in Example [T] with the general method. All the 
splines we use with the general method are B-splines. For this example, the knots are from 
5 to 48 at step 0.05, and we put one price level at the center of each knot interval. 

We solved fl4.3p with IPOPT using the full monotonicity constraint, and obtained the optimal 
value p = 0.0048, which is sufficiently close to 0. Figure (521 shows the spline approximation 
of the SFE. We see that when the mesh is hne enough, splines are quite capable at handling 
nonsmoothness and even discontinuities of the functions. 


Comparing Figure 15721 with Figure l5.1b[ we see that both methods are able to hnd the equi- 
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(a) Solutions with different values of t (b) Approximated SFE 

Figure 5.1: Least squares approximation 



Figure 5.2; General method approximation 

librium for markets of asymmetric duopoly with constant marginal costs, and the solutions 
are both of high precision. However, no matter in terms of the computational time, or in 
terms of the tools, solving a least squares problem is far easier than solving a highly nonlin¬ 
ear large scale optimization problem. Thus for this type of problems, the specialized least 
squares method is certainly more preferable. 

Example 3. In this example we compare the effects of the two monotonicity constraints, 
and the results with different mesh sizes. The example is taken from Anderson and Hu 
|16) . which has three firms with cost functions Ci{q) = 5q + 0.8q‘^, C 2 {q) = 8q + 1.2g^ and 
Csi^q) = 12q + 2.3q‘^, and capacities Capi = 11, Cap 2 = 8 and Gaps = 55, respectively. The 
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(a) pointwise monotonicity constraint 


(b) full monotonicity constraint 


Figure 5.3: Comparison of monotonicity constraints 


demand function is D{p) = — 0.5p + e. 


First we compare the pointwise monotonicity with the full monotonicity constraints. The 
knots we use are from 5 to 54 at step 0.5, and we put a price level at the center of each 
knot interval. We solve the problem with CONOPT; The optimal value of p is 1.6 x 
if we use the pointwise constraint, and it is 0.002 if we use the full constraint. Although 
the pointwise constraint gives a smaller p, it does not necessarily mean that it is the better 
choice. Figure 15.31 is a comparison of the results at the low price level, where images are 
magnihed. Theoretically, we know that when 5 < p < 8, Firm 1 is the monopoly, thus 
'Si(p) = “ 5), and we also know that si should have a jump at p = 8. We see that 

the solution given by the full constraint fj5.3bD is closer to the true Si than the solution 
given by the pointwise constraint f|5.3ap is. Therefore, despite a larger value of p, the full 
monotonicity constraint is in fact more robust. 

We also see from Figure 15.31 that although the full monotonicity constraint gives a better 
approximation, it is still not close enough to the true equilibrium. This is due to the hneness 
of the mesh, and as we make the mesh hner, the result will be better. As an illustration, 
we reduce the knot interval from 0.5 to 0.1, and again put one price level at the center 
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(a) |r| = 0.5 (b) |t| = 0.1 


Figure 5.4: Comparison of mesh sizes 



Figure 5.5: Approximation with |r| =0.1 


of each new knot interval. Keep the full monotonicity constraint and solve the problem, 
CONOPT gives a new result with p = 0.00017. Figure [531 compares the result of the hner 
approximation fj5.4bp with the previous coarser approximation fl5.4ap . It is apparent that 
the precision has signihcantly improved. Figure l575l is the full plot of the hner approximation 
for 5 < p < 54. 


For an even better approximation, one can always make the knot intervals smaller. However, 
as the number of knots increases, the problem will eventually become too large to solve. 
One way to improve the precision while keeping the problem size tractable is to use the 
information we obtained from a coarser approximation. From Figure 15751 we can see that 
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for 15 < p < 40, the supply functions are very smooth with little fluctuation, which means 
that a few pieces of quadratic polynomials are good enough to approximate them. Therefore 
we can, for example, set one knot at every 0.05 unit for 5 < p < 15, and at every 5 units 
for 15 < p < 40, which allows us to save hundreds of knots that would incur thousands of 
constraints. 


6 Conclusions 

One of the reasons why finding SFE has been so difficult is that the supply functions do 
not have specific forms, so all the non-decreasing functions (bounded by capacity) have to 
be considered. To find these free-form functions, parameterization is almost inevitable, and 
splines, due to their flexibility, are arguably the best way of parameterization for the purpose 
of approximation. For duopolies with constant marginal costs, we found that the first order 
conditions are linear in the spline coefficients, allowing us to approximate the solutions of the 
ODE system by solving a least squares problem. We proved that when the demand can be 
sufficiently low with a positive probability, the solution space of the least squares problem is 
exactly the solution space of the ODE system. And since least squares problems are so easy 
to solve, we can obtain solutions of high precision by using very fine mesh, while still solving 
it fast. The solutions have a clean form, which allows us to find the equilibrium easily by 
searching for the supply functions that reach the capacities smoothly. 

We also used splines to improve the general purpose method given in Anderson and Hu 
|16) . Both their original method and our proposal should be equally accurate, but the use 
of splines enabled us to significantly reduce the number of decision variables and constraints 
used in the auxiliary NLP, thus making the problem easier to solve, and without the need 
of human intervention. 
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The solutions of both the specialized and the general purpose methods are proved to converge 
uniformly to the SFE. We also provided numerical examples to demonstrate the use of these 
methods, and the solutions are precise and reflect the theoretical properties of SFEs that we 
developed throughout the paper. 
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